Los modelos lineales son ampliamente utilizados en acuicultura para explicar, modelar o predecir la relación lineal de una variable respuesta \(Y\) con una o mĆ”s variables predictoras \(X_1,ā¦,X_p\).
\[Y_{i} = \beta_{0} + \beta_{1} X_{i1} + \beta_{2} \beta_{i2} + ... + \beta_{p} X_{ip} + \epsilon_{i}\]
Supuestos: \(Y_i\) se distribuye normalmente
Si p = 1, el modelo es una regresión lineal simple.
Si p > 1, el modelo es una regresión lineal múltiple.
Si p > 1 y alguna variable predictora es categórica, el modelo se denomina ANCOVA.
Prueba de hipótesis del coeficiente de regresión y el intercepto
La hipótesis nula en ambos casos es que el coeficiente del intercepto (\(β_0\)) y de la pendiente (\(β_1\)) son iguales a 0.
\(H_0:β_0 = 0\) y \(H_0:β_1 = 0\)
El Coeficiente de determinación (R²) se calcula como el cuadrado del coeficiente de correlación de pearson. Este nos indica que tan buena es la predicción que el modelo hace de los datos.
El R cuadrado ajustado (R² ajustado) nos dice qué porcentaje de la variación de la variable dependiente es explicado por la o las variables independientes de manera conjunta. En un modelo de regresión lineal, el R² ajustado es una medida de bondad que considera el número de variables existentes en el modelo.
\[R^2_{ajust} =1-(1-R^2)\frac{n-1}{n-p-1}\]
donde:
\(n\) = tamaƱo de la muestra
\(p\) = cantidad de variables predictoras en el modelo
Prueba de hipótesis del modelo completo
La hipótesis nula es si los coeficientes son iguales a 0.
\(H_0:β_j = 0\) ; \(j = 1, 2,...,k\)
Los objetivos de aprendizaje de esta guĆa son:
1. - Realizar anƔlisis de datos usando modelos lineales.
2. - Realizar grƔficas avanzadas con ggplot2.
3. - Elaborar un reporte dinƔmico en formato html con Rmarkdown.
En vertebrados, la miostatina (MSTN) actua como es un represor del crecimiento muscular. Muchos trabajos han descrito que mutacioones en ese gen producen un fenotipo llamado de doble músculo caracterizado por un crecimiento exacerbado del músculo esquelético. En este ejercicio trabajaremos con un set de datos simulado de perros de carrera que poseen este distintivo fenotipo Adaptado de Mosher et al. 2007
Figura 1: Comparación de animales con diferentes genotipos del gen Miostatina: (A) +/+; (B) 1 mutante con alelo cys ā stop (mh/+); (C) Dos copias de la mutación cys ā stop (mh/mh).
La variable respuesta se denomina mass-to-height ratio y se expresa en Kg/cm. Para este set de datos existen dos variables predictoras que se deberƔn analizar: 1) Geno: que corresponde al genotipo de los animales; 2) Sex: que corresponde al sexo.
La variable Geno tiene 3 niveles denominados 0, 1 y 2 equivalentes a los genotipos +/+, +/mh y mh/mh donde 0, 1 y 2 representan el número de alelos mh. Esta forma numérica de representar los genotipos es conveniente para analizar la variable Geno como una variable cuantitativa discreta 0, 1 y 2 o como un factor (ej. factor(Geno)) según sea nuestro interés. La variable Sex tiene dos niveles male y female.
Elabore un documento .Rmd y configure su reporte para exportar en .html. Instale los paquetes readxl, ggplot2, dplyr y multcomp para el anƔlisis de los datos.
knitr::opts_chunk$set(echo = TRUE)
library(readxl)
library(ggplot2)
library(dplyr)
library(multcomp)
Ejecute cada uno de los siguientes ejercicios en uno o mÔs bloques de códigos diferentes. Sea ordenado y documente su reporte adecuadamente.
Importe el set de datos myostatin.deletion.data.xlsx y transforme solo la variable Sex a factor, luego realice un anƔlisis exploratorio de datos.
Incluya:
a). Resumen estadĆstico de todas las variables.
snp.data <- read_excel("myostatin.deletion.data.xlsx")
snp.data$Sex <- as.factor(snp.data$Sex)
summary(snp.data)
## animal Geno Sex mass-to-height ratio
## Min. : 1 Min. :0.0000 female:494 Min. :0.2009
## 1st Qu.:174 1st Qu.:0.0000 male :199 1st Qu.:0.2919
## Median :347 Median :1.0000 Median :0.3401
## Mean :347 Mean :0.9769 Mean :0.3590
## 3rd Qu.:520 3rd Qu.:2.0000 3rd Qu.:0.4100
## Max. :693 Max. :2.0000 Max. :0.5993
b). Tabla de frecuencia combinada Geno por Sex, esto es importante para ver si los datos son o no balanceados.
table(snp.data$Sex, snp.data$Geno)
##
## 0 1 2
## female 155 80 259
## male 120 79 0
tapply(snp.data$`mass-to-height ratio`, factor(snp.data$Geno), mean)
## 0 1 2
## 0.2829893 0.3263597 0.4597127
tapply(snp.data$`mass-to-height ratio`, factor(snp.data$Sex), mean)
## female male
## 0.3756034 0.3177426
c). Histograma de la variable respuesta mass-to-height ratio.
ggplot(snp.data, aes(x=`mass-to-height ratio`))+
geom_histogram(color="darkblue", fill="lightblue", bins = 12)
a). Realice una grĆ”fica de dispersión (geom_point) de weight en función de Geno, incluya el comando geom_smooth(method=lm) para agregar la lĆnea de regresión a la grĆ”fica.
p <- ggplot(snp.data, aes(x = Geno, y = `mass-to-height ratio`))
p + geom_point() + xlab("Reference allele count") + geom_smooth(method=lm)
b). Realice un anÔlisis de regresión lineal para investigar la asociación entre weight y Geno usando las funciones lm(), summary(). Interprete los resultados del modelo lineal y responda las siguientes preguntas. ¿Qué representa el estimador de Geno y del intercepto?. ¿La pendiente es significativamente distinta de cero?. ¿Y el intercepto?
lm.geno <- lm(`mass-to-height ratio` ~ Geno, data = snp.data)
summary(lm.geno)
##
## Call:
## lm(formula = `mass-to-height ratio` ~ Geno, data = snp.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.110802 -0.044236 -0.006808 0.044296 0.150224
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.272969 0.003434 79.49 <2e-16 ***
## Geno 0.088052 0.002615 33.67 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06041 on 691 degrees of freedom
## Multiple R-squared: 0.6213, Adjusted R-squared: 0.6208
## F-statistic: 1134 on 1 and 691 DF, p-value: < 2.2e-16
c). Realice la predicción de la media poblacional de weight para los genotipos AA, AB y BB y la predicción del peso weight de tres nuevos individuos con genotipo AA, AB y BB. ¿Qué tan buena son las predicciones obtenidas? Analice el R² ajustado y los intervalos de confianza.
predict.lm(lm.geno, newdata=data.frame(Geno=c(0,1,2)), interval="confidence")
## fit lwr upr
## 1 0.2729690 0.2662268 0.2797112
## 2 0.3610212 0.3565142 0.3655282
## 3 0.4490734 0.4421530 0.4559937
predict.lm(lm.geno, newdata=data.frame(Geno=c(0,1,2)), interval="prediction")
## fit lwr upr
## 1 0.2729690 0.1541725 0.3917655
## 2 0.3610212 0.2423306 0.4797118
## 3 0.4490734 0.3302666 0.5678801
d). Realice una anova del modelo lineal con la función anova() (supone varianzas homogeneas) y compare con una anova que no supone varianzas homogeneas oneway.test()
# anova suponiendo varianzas iguales
anova(lm.geno)
## Analysis of Variance Table
##
## Response: mass-to-height ratio
## Df Sum Sq Mean Sq F value Pr(>F)
## Geno 1 4.1373 4.1373 1133.8 < 2.2e-16 ***
## Residuals 691 2.5215 0.0036
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# anova suponiendo varianzas diferentes
oneway.test(`mass-to-height ratio` ~ Geno, data = snp.data)
##
## One-way analysis of means (not assuming equal variances)
##
## data: `mass-to-height ratio` and Geno
## F = 520.26, num df = 2.00, denom df = 427.93, p-value < 2.2e-16
e) El siguiente modelo lineal permite realizar pruebas de contraste entre los tres diferentes genotipos.
lm(weight ~ -1 + factor(Geno), data = snp.data)
Compare los Betas y el \[R^2\] de este modelo con el realizado en el ejercicio b.
lm.geno2 <- lm(`mass-to-height ratio` ~ -1 + factor(Geno), data = snp.data)
summary(lm.geno2)
##
## Call:
## lm(formula = `mass-to-height ratio` ~ -1 + factor(Geno), data = snp.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.119005 -0.044944 -0.000312 0.040813 0.139584
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## factor(Geno)0 0.282989 0.003461 81.75 <2e-16 ***
## factor(Geno)1 0.326360 0.004552 71.69 <2e-16 ***
## factor(Geno)2 0.459713 0.003567 128.89 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0574 on 690 degrees of freedom
## Multiple R-squared: 0.9763, Adjusted R-squared: 0.9762
## F-statistic: 9478 on 3 and 690 DF, p-value: < 2.2e-16
summary(lm.geno)
##
## Call:
## lm(formula = `mass-to-height ratio` ~ Geno, data = snp.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.110802 -0.044236 -0.006808 0.044296 0.150224
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.272969 0.003434 79.49 <2e-16 ***
## Geno 0.088052 0.002615 33.67 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06041 on 691 degrees of freedom
## Multiple R-squared: 0.6213, Adjusted R-squared: 0.6208
## F-statistic: 1134 on 1 and 691 DF, p-value: < 2.2e-16
f) Investigue y aplique los comandos contrMat() y glht para realizar una prueba de contrastes con ajuste de bonferroni para comparaciones mĆŗltiples.
# Elabora matriz de contrastes para el factor Geno
contrastes = contrMat(table(snp.data$Geno), type="Tukey")
# Realiza comparaciones multiples
mc2 = glht(lm.geno2, linfct =contrastes)
summary(mc2, test=adjusted("bonferroni"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = `mass-to-height ratio` ~ -1 + factor(Geno), data = snp.data)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 - 0 == 0 0.043370 0.005719 7.584 3.26e-13 ***
## 2 - 0 == 0 0.176723 0.004970 35.556 < 2e-16 ***
## 2 - 1 == 0 0.133353 0.005783 23.059 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)
a). Realice una grĆ”fica de dispersión (geom_point) de mass-to-height ratio en función de Geno y Sex, incluya el comando geom_smooth(method=lm) para agregar la lĆnea de regresión a la grĆ”fica.
q <- ggplot(snp.data, aes(x = Geno, y = `mass-to-height ratio`, shape=Sex, color=Sex))
q + geom_point() + xlab("Reference allele count") + geom_smooth(method=lm)
b). Realice un anÔlisis de covarianza para investigar la asociación entre mass-to-height ratio, Geno y Sex usando las funciones lm(), summary().
lm.geno.sex <- lm(`mass-to-height ratio` ~ Geno + Sex, data = snp.data)
summary(lm.geno.sex)
##
## Call:
## lm(formula = `mass-to-height ratio` ~ Geno + Sex, data = snp.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.107264 -0.046971 -0.006527 0.041388 0.151325
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.264638 0.004396 60.205 < 2e-16 ***
## Geno 0.091667 0.002864 32.006 < 2e-16 ***
## Sexmale 0.016714 0.005555 3.009 0.00272 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06006 on 690 degrees of freedom
## Multiple R-squared: 0.6262, Adjusted R-squared: 0.6251
## F-statistic: 578 on 2 and 690 DF, p-value: < 2.2e-16
c). Interprete los resultados del modelo lineal.
a). Realice grÔficas del tamaño de los efectos plot.design() y de interacción interaction.plot de Geno y Sex sobre mass-to-height ratio.
plot.design(snp.data$`mass-to-height ratio` ~ factor(snp.data$Geno) + factor(snp.data$Sex))
interaction.plot(factor(snp.data$Geno), factor(snp.data$Sex), snp.data$`mass-to-height ratio`, mean)
b). Realice un anÔlisis de varianza con interacción para investigar la asociación entre weight, Geno y Sex usando las funciones lm(), summary().
lm.qtl <- lm(`mass-to-height ratio` ~ factor(Geno) * factor(Sex), data = snp.data)
summary(lm.qtl)
##
## Call:
## lm(formula = `mass-to-height ratio` ~ factor(Geno) * factor(Sex),
## data = snp.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.119005 -0.041228 0.000675 0.040398 0.139584
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.271157 0.004495 60.321 < 2e-16 ***
## factor(Geno)1 0.034509 0.007704 4.479 8.78e-06 ***
## factor(Geno)2 0.188556 0.005683 33.177 < 2e-16 ***
## factor(Sex)male 0.027117 0.006805 3.985 7.47e-05 ***
## factor(Geno)1:factor(Sex)male 0.014534 0.011185 1.299 0.194
## factor(Geno)2:factor(Sex)male NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05597 on 688 degrees of freedom
## Multiple R-squared: 0.6764, Adjusted R-squared: 0.6745
## F-statistic: 359.5 on 4 and 688 DF, p-value: < 2.2e-16
c). Interprete los resultados del modelo lineal.
a). Evalue mormalidad del modelo con interacción del ejercicio 5 usando métodos grÔficos y prueba de shapiro.
shapiro.test (residuals (lm.qtl))
##
## Shapiro-Wilk normality test
##
## data: residuals(lm.qtl)
## W = 0.99161, p-value = 0.0005773
plot(lm.qtl, which = 2)
b). Evalue homogeneidad de varianzas usando el comando plot.
plot(lm.qtl, which = 1)